# -*- coding: utf-8 -*-
"""
Created on Thu Sep 27 12:07:21 2012
This module reads the engine database and processes the data to solve for
engine performance characteristics

The module provides the following functions
    - database i/o
    - temperature correction
    - partial throttle analysis
    - power, fuel flow, and RPM at input flight conditions and throttle

@author: daniel
"""
import dbTools
from scipy.interpolate import interp1d, Rbf
import numpy as np
from FlightConditions import ISAtmosphere as atm
import paths
import matplotlib.pyplot as plt


def load(name):
    db          =paths.Database()
    engineDB    =dbTools.loadDB(db.engine)
    engineSheet =engineDB.selectByName("ROTAX_912ULS")
    eng         =pistonEngine(engineSheet,False)
    return eng

class pistonEngine(object):
    """    
    Contains functions for reading the engine database and generating
    spline curves, and returning performance data
    
    Parameters
    ----------
    
    entryName : string
        name of engine in the database to be loaded
    
    
    Examples
    --------
    
    >>> import paths
    >>> db          =paths.Database()
    >>> engineDB    =dbTools.loadDB(db.engine)
    >>> engineSheet =engineDB.selectByName("ROTAX_912ULS")
    >>> eng         =pistonEngine(engineSheet,False)
    >>> print eng.get_availablePower(1.05,2000.,100)
    26823.7503111
    >>> print eng.sfc[0]
    4.38107638889e-08
    >>> print eng.get_idle_parameters()
    (32000.0, 0.0014019444444444444)
    """
    def __init__(self,engineInputSheet,suppressWarnings=True):
        self._inputSheet  =engineInputSheet
        self.entryName    =str(self._inputSheet.name)
        self.name         =str(self._inputSheet.cell(0,1).value)
        self._readData()
        self._create_curves()
        self.suppressWarnings=suppressWarnings
    def _readData(self):
        """
        reads engine database excel file.
        
        Parameters
        ----------
        
        gearboxRatio  : float
            engine gearbox ratio
        mass : float, kg
            engine mass
        length : float, m
            engine length
        width : float, m
            engine width
        height : float, m
            engine height
        flatRateAlt : float, m
            engine flat rate altitude
            The altitude at which engine performance begins to decrease due to
            reducing atmospheric density
        RPM : array
            array of rpm from db
        torque : array
            array of torque at rpm
        fuelFlow : array
            array of fuel flow at rpm
        fuelDensity : float
            density of fuel for conversion kg/h, l/h
        """
        db=dbTools.readDB(self._inputSheet)
        self.gearboxRatio =db.readRow(1,1)
        self.mass         =db.readRow(2,1)
        self.length       =db.readRow(3,1)
        self.width        =db.readRow(4,1)
        self.height       =db.readRow(5,1)
        self.flatRateAlt  =db.readRow(6,1)
        self.RPM          =db.readRow(7,1)
        self.power        =db.readRow(8,1)*1000
        self.torque       =db.readRow(9,1)
        self.fuelFlow     =db.readRow(10,1)
        self.fuelDensity  =db.readRow(11,1)
        self.sfc          =self.fuelFlow*self.fuelDensity/3600.0/self.power
        
    def _create_curves(self):
        """
        creates interpolation curves using cubic spline interpolation of tabulated data.
        
        Parameters
        ----------
        
        rpmPowerCurve : spline
            curve used to determine rpm at given power setting
        powerCurve : spline
            rpm vs. power spline
        torqueCurve : spline
            rpm vs. torque spline
        fuelCurve : spline
            rpm vs. fuelFlow spline
        """
        self.rpmPowerCurve = interp1d(self.power,self.RPM,'cubic')
        self.powerCurve    = interp1d(self.RPM,self.power,'cubic')
        self.torqueCurve   = interp1d(self.RPM,self.torque,'cubic')
        self.fuelCurve     = interp1d(self.RPM,self.fuelFlow,'cubic')
        self.sfcCurve      = interp1d(self.RPM,self.sfc,'cubic')
        
    def density_correction(self,density):
        """
        Density correction considering flat rate altitude.
        """
        c            =0.132  
        flatRateAtm  =atm(self.flatRateAlt)
        rho0         =flatRateAtm.density
        if rho0 <= density:
            densCorr =1.0
        else:
            densCorr =(1.0+c)*(density)/rho0-c
        return densCorr
    def get_idle_parameters(self):
        """
        Returns idle power and fuel flow
        """
        power = self.power[0]
        fuelFlow = self.fuelFlow[0]*self.fuelDensity / 3600.0
        return power,fuelFlow

    def get_intakeDensity(self,density,powerSetting):
        return density*powerSetting/100.0
    
    def get_availablePowerCurve(self,density,powerSetting):
        powerCurve=self.power*self.density_correction(self.get_intakeDensity(density,powerSetting))
        return interp1d(self.RPM,powerCurve)

    def get_availablePower(self,density,RPM,powerSetting):
        """
        Returns available power in Watt
        
        Parameters
        ----------
        density : float, kg/m**3
            air density
        RPM : float
            engine rpm
        powerSetting : float [0;100]
            engine power setting
        """
        RPM=self._checkRPM(RPM)
        powerSetting=self._checkPowerSetting(powerSetting)
        powerCurve=self.get_availablePowerCurve(density,powerSetting)
        return powerCurve(RPM)

    def get_fuelFlow(self,density,RPM,powerSetting):
        """
        Returns fuel flow in kg/sec
        
        Parameters
        ----------
        density : float, kg/m**3
            air density
        RPM : float
            engine rpm
        powerSetting : float [0;100]
            engine power setting
        """
        RPM=self._checkRPM(RPM)
        powerSetting=self._checkPowerSetting(powerSetting)
        sfc=self.get_sfc(RPM)
        Pa =self.get_availablePower(density,RPM,powerSetting)
        fuelFlow=sfc*Pa        
        return fuelFlow

    def get_sfc(self,RPM):
        """ 
        returns SFT at given rpm
        """
        RPM=self._checkRPM(RPM)
        return self.sfcCurve(RPM)

    def get_maxRPM(self):
        """ 
        returns maximum rpm of the engine
        """
        return np.max(self.RPM)
        
    def get_minRPM(self):
        """ 
        returns minimum rpm of the engine
        """
        return np.min(self.RPM)
        
    def _checkRPM(self,RPM):
        minRPM=self.get_minRPM()
        maxRPM=self.get_maxRPM()
        if RPM < minRPM:RPM=minRPM
        if RPM > maxRPM:RPM=maxRPM
        return RPM
        
    def _checkPowerSetting(self,powerSetting):
        minPWR=10.
        maxPWR=100.
        if powerSetting < minPWR:
            powerSetting=minPWR
        if powerSetting > maxPWR:
            powerSetting=maxPWR
        return powerSetting


class EngineExperimental:
    def __init__(self,name):
        self.pth = paths.Database()
        self.atmSL = atm(0.0)
        self._read_xls(name)
    def _read_xls(self,name=None,dbName=None):
        if dbName==None:
            dbName = self.pth.engineExp
        if name==None:
            name = 'ROTAX_912'
        engDb = dbTools.loadDB(dbName)
        engSh = engDb.selectByName(name)
        db = dbTools.readDB(engSh)
        idx = db.findHeader('SETTING')
        self.rpmLb = np.zeros(len(idx))
        self.rpmUb = np.zeros(len(idx))
        psetList   = np.zeros(len(idx))
        for i,j in enumerate(idx):
            rpm = db.readRow(j+1,1,True)
            psetList[i] = db.readRow(j,1,False)
            pset = np.ones(len(rpm))*psetList[i]
            self.rpmLb[i] = rpm[0]
            self.rpmUb[i] = rpm[-1]
            if i==0:
                self.powerSetting = pset
                self.rpm = rpm
                self.power = db.readRow(j+2,1,True)
                self.fuelFlow = db.readRow(j+3,1,True)
            else:
                self.powerSetting = np.hstack([self.powerSetting,pset])
                self.rpm = np.hstack([self.rpm,rpm])
                self.power = np.hstack([self.power,db.readRow(j+2,1,True)])
                self.fuelFlow = np.hstack([self.fuelFlow,db.readRow(j+3,1,True)])
            self.rps = self.rpm / 60.0
        self.power = self.power * 745.699872 # HP to Watt
        self.powerSurf = Rbf(self.powerSetting,self.rpm,self.power)
        self.fuelFlowSurf = Rbf(self.powerSetting,self.rpm,self.fuelFlow)
        self.rpmLbCurve = interp1d(psetList,self.rpmLb,'cubic')
        self.rpmUbCurve = interp1d(psetList,self.rpmUb,'cubic')
        self.powerSettingList = psetList
    def _check_pset(self,pset):
        pMax = max(self.powerSettingList)
        pMin = min(self.powerSettingList)
        if pset>pMax:
            pset = pMax
        elif pset<pMin:
            pset = pMin
        return pset
    def _check_rpm(self,rpm,pset):
        rpmMax = self.get_rpm_max(pset)
        rpmMin = self.get_rpm_min(pset)
        if rpm>rpmMax:
            rpm = rpmMax
        elif rpm<rpmMin:
            rpm = rpmMin
        return rpm
    def get_rpm_min(self,pset):
        pset = self._check_pset(pset)
        return self.rpmLbCurve(pset)
    def get_rpm_max(self,pset):
        pset = self._check_pset(pset)
        return self.rpmUbCurve(pset)
    def _density_correction(self,density):
        return density / self.atmSL.density
    def get_availablePower(self,rpm,pset,density=None):
        if density==None:
            corr = 1.0
        else:
            corr = self._density_correction(density)
        pset = self._check_pset(pset)
        rpm = self._check_rpm(rpm,pset)
        return corr*self.powerSurf(pset,rpm)
    def _corr_function(self,x,lb,ub,w=1.0):
        p1 = w*(-lb+x)
        p2 = w*(ub-x)
        threshold = 500.0
        if p1>threshold:
            yl = 0.0
        elif p1<-threshold:
            yl = -1.0
        else:
            yl = -1.0 / (1.0+np.exp(p1))
        if p2>threshold:
            yu = 0.0
        elif p2<-threshold:
            yu = -1.0
        else:
            yu = -1.0 / (1.0+np.exp(p2))
        return yl+yu
    def get_availablePower_constr(self,rpm,pset,density=None):
        rpmLim = [self.get_rpm_min(pset),self.get_rpm_max(pset)]
        Pav = self.get_availablePower(rpm,pset,density)
        corr = self._corr_function(rpm,rpmLim[0],rpmLim[1],0.05)
        return Pav + corr*Pav
    def get_fuel_flow(self,rpm,pset):
        pset = self._check_pset(pset)
        rpm = self._check_rpm(rpm,pset)
        return self.fuelFlowSurf(pset,rpm)
        
def runTest():
    import paths
    db          =paths.Database()
    engineDB    =dbTools.loadDB(db.engine)
    engineSheet =engineDB.selectByName("ROTAX_912ULS")
    eng         =pistonEngine(engineSheet,False)
    print eng.get_availablePower(1.05,2000.,100)
    print eng.sfc[0]
    print eng.get_idle_parameters()

def runTest2():
    from numpy import linspace,zeros, array
    
    rpmExp1 = array([1700,2000,2300,2600])
    QExp1 = array([5.38,5.26,4.39,3.83])
    Pexp1 = rpmExp1*QExp1
    rpmExp = array([2300,2600,2900,3200,3500])
    QExp = array([7.89,7.74,7.46,7.24,6.69])
    Pexp = QExp*rpmExp
    rpmExp2 = array([4400,4700,5000,5300,5600,5800])
    QExp2 = array([12.28,12.35,12.25,12.05,11.76,11.50])
    Pexp2 = rpmExp2*QExp2
    
    eng = load('ROTAX_912ULS')
    n = 15
    #pSet = linspace(10.0,100.0,n)
    pSet = [20.,50.0,75.0,95,100.0]
    rpm = linspace(eng.get_minRPM(),eng.get_maxRPM(),n)
    crs = atm(1500.0)
    plt.figure(1)
    plt.grid(True)
    plt.hold(True)
    leg = list()
    for p in pSet:
        pCurve = eng.get_availablePowerCurve(crs.density,p)
        plt.plot(rpm/2.43,pCurve(rpm)/1e3,'o-')
        leg.append('%.0f%% throttle'%p)
#    plt.figure(2)
#    plt.hold(True)
#    plt.plot(rpmExp,Pexp)
#    plt.plot(rpmExp1,Pexp1)
#    plt.plot(rpmExp2,Pexp2)
    plt.legend(leg,'upper left')
    plt.xlabel('rpm')
    plt.ylabel('power, kW')
    plt.show()

def run_test3():
    name = 'ROTAX_912'
    expEngine = EngineExperimental(name)
    pset = 75.
    rpm = np.linspace(0.0,5000.0,1000)
    p = np.zeros(len(rpm))
    for i,r in enumerate(rpm):
        p[i] = expEngine.get_availablePower_constr(r,pset,1.05)
    plt.figure()
    plt.plot(rpm,p)
    plt.show()

if __name__=="__main__":
    runTest()